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ABSTRACT 


The method for combining geometrical data with satellite dy- 
namical and gravimetry data for the solution of geopotential and 
station location parameters is discussed. Geometrical tracking 
data (simultaneous events) from the global network of BC-4 
stations are currently being processed in a solution that will 
greatly enhance the geodetic world system of stations. Previ- 
ously the stations in Goddard Earth Models have been derived 
only from dynamical tracking data. In this paper a linear re- 
gression model is formulated for combining the data, based upon 
the statistical technique of weighted least squares. Reduced 
normal equations, independent of satellite and instrumental pa- 
rameters, are derived for the solution of the geodetic parame- 
ters. Exterior standards for the evaluation of the solution and 
for the scale of the earth's figure are discussed. 
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GEOMETRICAL GEODESY TECHNIQUES 
IN GODDARD EARTH MODELS 


I. INTRODUCTION 

A matrix model employing a reduced form of the general least squares adjust- 
ment process is developed. The model provides a method for combining geo- 
metrical BC-4 optical data with satellite dynamic and gravimetric data into a 
general geodetic solution. The geodetic solution consists of a geocentric system 
of station coordinates and spherical harmonic coefficients for the geopotential. 

A previous solution, Lerch et al (1972), for a Goddard Earth Model (GEM 4) con- 
tained 514 geodetic parameters but did not include any geometric data. 

A map of station locations is presented in Figure 1, illustrating the distribution 
of 45 BC-4 stations associated with the geometric data and 61 stations associated 
with dynamic data from electronic, laser, and optical tracking systems used pre- 
viously in GEM 4. Figure 2 illustrates local datum ties between the dynamic and 
geometric stations and BC-4 baselines obtained from the geodimeter and tellu- 
rometer systems that are to be employed in the new combination solution. 


n. BASIC DATA SYSTEMS AND REPRESENTATION 
1. Data Systems 

Four basic geodetic data systems are employed in the combination solution and 
are listed with their associated solution parameters in Table 1. _These systems_ 
are the BC-4 geometric (B), gravimetry (G), dynamic satellite (D), and survey (S). 

In the survey data system S, the baselines and ties were illustrated in Figure 2. 
The ties correspond to relative position coordinates for nearby geometrical and 
dynamical stations and are treated as observations with statistical errors based 
upon local survey accuracy. Similarly the baseline distances are treated as 
statistical observations. 

The satellite related parameters associated with the BC-4 system B and the 
satellite dynamic system D require preliminary processing for initial estimates 
of the parameters. The initial processing will be discussed in a later section. 
There are some 20,000 satellite position parameters associated with 1100 BC-4 
events of two, three, and four stations observing the satellite simultaneously at 
6 to 7 reduced points of satellite position for each event. There are some 350 
weekly arcs of satellite dynamic data on 27 satellites, where some 400,000 ob- 
servations of electronic, laser, and optical tracking data have been processed. 
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Figure 1. Station 'Locations 




Figure 2. Baselines and Station Ties (Survey Data) 




Table 1 

Basic Data Systems and Solution Parameters 


Symbol 

Data System 

Solution Parameters 
(unknowns) 

Unknown 

Symbol 

5 

BC-4 geometric 

station location coordinates 

X 



satellite position components 
for each geometrical point j 


G 

Gravimetry 

potential coefficients 
(spherical harmonic) 

c 

D 

Dynamic satellite 

potential coefficients 

c 



station coordinates 
satellite orbital elements 

z 



and tracking system parameters 
associated with each satellite 
arc i 

Pi 

S 

Survey 

station coordinates connecting 



(Baselines and ties) 

BC-4 baselines x b in x 

X 



station coordinates connecting 
ties x t in x, z t in z 

X, z 


The satellite dynamic parameters consist of six orbital elements and modeled 
force parameters at a given epoch on each weekly arc. Tracking system parame- 
ters for certain electronic systems have been included in the modeling and are 
associated with the initial processing on a weekly arc. The gravimetry data 
system G consists of a global distribution of 5° equal area blocks of mean gravity 
anomalies, Rapp, (1972). 

Because of the large number of satellite parameters in the systems B and D, a 
least squares matrix model is developed that will reduce the matrix to a form 
containing just the geodetic parameters. The reduced form will then be suitable 
for computational solution. 


2. Representation of Data Systems 

A total data system C, which will be used to encompass the four basic data 
systems, is represented as 
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C : (0, C, v; y) 


( 1 ) 


where each symbol denotes a column vector as follows: 

O — observations 

C — computed quantities corresponding to the observations 

v — observation residuals (v = O - C) 

y — the solution parameters or unknowns (see Table 1) 

The column vectors are further defined in Section III where the method for the 
solution is developed. The four basic data systems, subsystems of C, are simi- 
larly represented and defined as in the above form (1), namely 


B : 

(0 B , B, b; 

X, q) 

(2) 

G : 

: (0 0 , G, g; 

c) 

(3) 

D : 

(Ob, D, d; 

C, z, p) 

(4) 

S : 

(O s , S, s; 

z, x) 

(5) 


where the solution parameters c, z, x, p = [p.] , q = £ qj ] are defined in Table 1 
for each of the data subsystems. Using the form above, data subsystems of D 
find b are represented for each dynamic satellite arc i and geometric satellite 
position point j, respectively as follows: 


D. 

1 

: (0 Di . D i , 

c, z, pj) 

(6) 

B. 

j 

: (0 Bj> B. , 

; x, q s ) 

(7) 
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IEL. DEVELOPMENT OF THE METHOD 

1. General Matrix Model for the Least Squares Solution 

A general matrix model of_the least squares adjustment process is developed for 
the complete data system C given in (1). The model will be subsequently employed 
to develop a reduced form of solution for the geodetic data subsystems. 

In the data system C the vector symbols, O, C, v, and the unknown solution y 
were defined under (1). The vector C = C (y) and the residual vector is 


v - 0 - C(y ) 


( 8 ) 


which in component form is 


[vj - [0„ ~ C n (y)] for n - 1 to N 


(9) 


where C n (y) is the computed quantity corresponding to the observation O n and is 
a function of the parameters in y. The linear condition equation for v in the 
least squares adjustment process is given by use of Taylor's expansion as 


v - v o - C y A y , 


( 10 ) 


where v = O - C (y o ), y = Ay + y c , y Q is a column vector of initial estimates for 
each component y k in y for k = 1 to K, and the matrix of partial (derivative) co- 
efficients 


C 


y 


~ BC n(y) 

. _ 


NXK 


( 11 ) 


in which the partial coefficient element lies in row n and column k and C n (y) is 
given under (9) above. C y is evaluated at y = y o . 

The least squares minimum condition is 


Q - v T W v v - minimum 


( 12 ) 
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where W v is a diagonal weight matrix with each element 


w_ 


for n - 1 to N 


( 13 ) 


in which a 2 is the error variance for the observation O . The minimum condi- 

n n 

tion for (12) is obtained from 


3Q _ 


f3 Ql 

= T 1 = 0 

L By kJ 


for k - 1 to K 


(14) 


which with use of (10) and (12) will give the normal equation 


C T W v = 0 

y v 


(15) 


or 


(C T W C )A y -C T W v =0. 

>■ y v y / J y v o 


(16) 


Solving (16) for Ay will give 


(C T W C ) _1 C T W v 

V y vy/ y v o 


(17) 


and the solution y is 


y - y 0 + A y- 


Under ideal conditions, where the modeling of C n (y) for O n is complete except 
for a random observation error which has normal distribution and variance crj 2 , 
the least squares solution for y may be shown to satisfy the maximum likelihood 
principle. See Anderson (1958). 



2. Reduced Form of the Normal Equations for the Geodetic Subsystems 


Using the column matrices defined for the total data system C under (1) and 
those for the geodetic data subsystems in (2) through (5), the following matrix 
partitioning is given for C in terms of its data subsystems: 


V 


V 


d 



B 


b 

c = 


II 

u 

1 

o 

II 

> 


°G 


G 


g 

_°s_ 


S 


s 


(18) 


Also partition 




c 


Ac" 



z 


A z 

\ 

y - 

X 

Ay = 

A x 

w 





g 


p 


Ap 

W s _ 

I 



I 

> 

1 


(19) 


where the solution parameters p = [p J , q = [q.] , x z, and c are defined in 
Table 1, the weight matrix W y was defined under (13) for C and W d , W b , W g , 

W s are defined as diagonal weight matrices respectively for the subsystem ob- 
servations as ordered in O above. When v and W v above are substituted into (12), 
Q then becomes 


Q = d T W d d + b T W b b + g T W g g + s T W s s (20) 

In terms of the variables associated with the data subsystem as given in Table 1, 
the linearized residual equations maybe expressed by Taylor's series for each 
subsystem as 
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d = d o " D c A c - D z A z - D p A p 

b = b o " B * A x - B q A q 

g = g c - G c A c 

s = s 0 - S x A x - S z A z 


( 21 ) 


where the subscript (o) corresponds to the initial estimate of the residual vector 
which is analogous to v o under (10), and the matrices of partial coefficients are 
defined analogous to C y under (11), for example, in D z above D would correspond 
to C and z to y. Differentiating Q with respect to each of the variables as ordered 
in (19) for y and setting the result equal to zero in order to obtain the minimum, 
the following normal equations will result with use of (21) for each of the vari- 
ables: 


D„ T W„ d + Gj W g g = 0 
BJ W d d + S, T W, s = 0 
BjW b b + SjW,s - 0 
D p tW d d = 0 
BjW b b - 0 


for c 
for z 

for x (22) 

for p 
for q 


Substituting the relations for d, b, g, and s in (21) into (22) will produce the 
normal matrix equations in terms of the unknown parameters and the initial 
residuals. This result between (21) and (22) is equivalent to the result between 
(15) and (16) for the total data system C. If C y is partitioned in terms of the 
partial coefficient matrices in (21) then by substituting its transpose, v from (18), 
and W v from (19) into (15) the normal equations (22) will result directly. 

The last two equations in (22) are used with d and b in (21) to express the satellite 
related parameters (Ap, A q) in terms of the geodetic parameters (Ac, Ax, Az). 
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When these results for Ap and Aqare substituted into d and b for the first three 
equations of (22), the reduced normal equations will then be obtained for the 
geodetic parameters. This process just described is carried out below. 

Proceeding in this manner then from (22) and d in (21) 

D p TW dd = D p T W d (d c - D c Ac - D z Az - D p A P ) = 0. (23) 


Solving (23) for Ap and substituting the result back into d will give 


d = P (d o - D c Ac - D z A z) (24) 

where 


P = I - D P 

P 

(25) 

P = (DjW d D)-*D;W d 

(26) 

Ap = P (d o - D c A c - D z A z) 

(27) 

Proceeding similarly as in (23) with 


B q T \b = 0 


then 


b = Q (b 0 - B x A x) 

(28) 

where Q and Aq may be derived as in (25) through (27). 


By use of d in (24), b in (28), and s and g in (21) the system (22) 
matrix form for the geodetic parameters as follows: 

is expressed in 
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Dec + Gcc 

Dcz 

0 


A c 


Deo 

+ Geo 

Dcz T 

Dzz + Szz 

S z x 


A z 

- 

Dzo 

+ Szo 

0 

Szx T 

Bxx + Sxx 

i 

_Ax_ 


Bxo 

+ Sxc_ 


( 29 ) 


where the above symbols are defined as follows: 


dynamic satellite system D 


Dec 


D T W D 

c c 


Dc z - 


D T W D 

c z 


Dz z - 


D t W D 

z z 


(30) 


Deo = D T Wd 


c o 


Dzo - 


d t w d 

z o 


w = 


W d p 


BC-4 geometric system B 


Bxx = 


B X T \ Q B x 


Bxo - B x t W b Q b o 


(31) 


gravimetric system G 


Gee = G c T W g G c 


Geo = G t W g 


C g 


(32) 


S z z - 


s T w s 

z s z 


survey system S 


S zx 


s 1 w s 

z s x 


Sxx 


S * T w s x 


(33) 


Szo = S T W s 


z s o 


Sxo - 


S T W s 

x s o 
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In the computations the results for D in (30) and B in (31) are expressed in terms 
of reduced results obtained from the observations on each satellite arc i with 
satellite parameters p A and on each satellite geometric pomt j with coordinates 
qj . This situation is represented by the data subsystems D^in (6) and Bj in (7). 
The data sets and associated partial coefficient matrices of D and B may be par- 
titioned in terms of these quantities for their respective subsystems to obtain the 
detailed results. 

The procedure for these results will be briefly derived. From Dj in (6) and D 
in (4) 

d i = °Di “ C i ( c > z > Pi) for each i = 1 to I, (34) 


d = Ed. ] D = [D. ] Ap = [Ap.] p = [p.] (35) 

Using Taylor's expansion for d. in (34) 

^ = d io -H ; Ac - ZAz-PiAp. (36) 

then, from d in (21) and since d = [dj , 

d o = 


D 


Defining the weight matrix W di for 0 Di and ordering it on the diagonal of W d for 
i = 1 to I and using the above results for (34) through (37), then (23) through (27) 
become for each arc i 
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0 


P i TW di d i = P i TW di( d io - HiAc - Z, Az - P.Ap.) = 
d i = P i( d io - ^ Ac - Z. Az) 

P i = I -P A ^ (38) 

P i = (PjWdiPi)’ 1 p i T w di 


A Pi = P i( d io - Hi A c - ZjAz) 
and the results in (30) become 


Dcc " J"! (HTW^j) 

i 

Dcz (HTW.Z.) 

i 

Dzz = £] (ZjWiZ.) 

Deo (HTW^) 

Dzo r L (H ^ idio) 

Wi =w dl p i 


i i 


where the sum ranges from i = 1 to I. 

Treating the system B. in (7) in a like manner to D. above, the result in (31) 
would become 


Bxx 


J 


L < X > T ", x i > 


Bxo - 


J 




(40) 


where X . corresponds to Z . , W. to W . , and b to d . 

J 1 } l jo i o 

The solution for the geodetic parameters is obtained through the inverse of the 
reduced normal matrix in (29). The inverse matrix also represents the variance- 
covariance matrix for the geodetic parameters, as the weighting for the observa- 
tions given in (13) is inversely equal to the variance of the observation errors. 
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3. Computational Procedure 


In practice the contributions to the reduced normal matrix for each of the geo- 
detic subsystems are computed separately and then combined as in (29). The 
terms for the reduced normal matrix for each of these four systems are given 
separately in equations (30) through (33)^ The_abovejpractme is similarly true 
in the computations for the subsystems of D and B. of B where their respec- 
tive matrix terms are identified and combined as in (39) and (40). 


The initial starting values for the geodetic parameters are obtained from a pre- 
vious geodetic solution, survey data, or through a separate process of data 
analysis. 

3.1 Satellite Dynamic System 

The initial observation residuals and dynamic satellite parameters for the system 
D, are obtained from a preliminary reduction of the observation data on a weekly 
satellite arc. In this process a bias parameter is modeled for certain electronic 
systems for the observation data on each satellite tracking pass. Since a large 
number of these parameters may occur within a weekly arc span of data, each 
bias parameter is eliminated through the back subsitution process at the end of 
the associated tracking pass. Numerical integration is employed for the satellite 
orbit and the preliminary reduction of the data on the weekly arc is carried out 
through the least square process of successive iterations. The normal matrix 
for the system is formed immediately after the preliminary reduction. 

3.2 Geometric Geodesy Technique 

In the computations for the BC-4 data system the reduced normal equations (31) 
may be obtained through the formulation of condition equations which are inde- 
pendent of satellite parameters. This technique is described in the appendix of 
the report, and it includes the constraint equations for the station coordinate 
ties and baseline distances from datum survey. In addition it provides for use 
of simultaneous MOTS and laser data, which has recently been analyzed by 
Reece and Marsh (1973) for stations in the area of the United States. 


IV. ANALYSIS OF SOLUTION AND GEODETIC RESULTS 

The geodetic solution may be tested and analyzed with the use of survey data in 
several areas of investigation. 

The mean sea level height (MSL) from station survey may be compared with the 
height of the station above the geoid as computed from the solution of the poten- 
tial coefficients and the geocentric station coordinates including the reference 


14 



ellipsoid, Lerch et al (1972). In view of the distribution of stations in Figure 1 
the dispersion of the differences between survey and computed values about a 
mean line may be analyzed. The MSL heights from survey are generally re- 
ported to be accurate to about a meter for these stations. Any significant offset 
in the mean line from zero differences may be associated with the scale for the 
reference ellipsoid, and thus the equatorial radius (a e ) of the reference ellipsoid 
may be adjusted. Any significant dispersion in the differences (including system- 
atic differences) may be analyzed in terms of geographical areas and in terms of 
various tracking systems such as the electronic, laser, and optical systems. 

The dispersion may be analyzed in terms of the geopotential model particularly 
in areas where the gravimetric measurements are not available. These analyses 
may be supported with separate tests for geoid heights and station coordinates 
as indicated below. 

Geoid heights computed from the potential model may be tested and analyzed in 
certain major sui-vey areas that contain astrogeodetic deflections of the vertical 
and detailed gravimetric data, Vincent et al (1972). 

Through adjustment for scale, orientation, and datum shift the station coordinates 
between datum survey and the solution may be analyzed. In such an analysis for 
a geocentric station solution by Marsh et al (1971), an rms agreement of 3.5 me- 
ters on station coordinate differences has been obtained for 20 stations on the 
North American Datum. 

With a solution derived from the geodetic data systems and use of the above 
analysis the following geodetic results may be obtained: 

1. A global geoid represented in terms of spherical harmonic coefficients. 

2. A world datum of station coordinates including an adjustment for scale, 
orientation, and datum shift for local datums. 

3. Mean equatorial radius(a e )for the Earth. 

4. Mean equatorial value of normal gravity (g e ) may be derived from the 
gravimetric data, Rapp (1972). 
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APPENDIX 


Reduced Normal Equations for Geometric Satellite Geodesy 
Including Constraints from Datum Survey 


The method developed here for the least squares normal equations is based 
upon the technique of formulating reduced condition equations where the satellite 
parameters have been eliminated. The data considered consists of simultaneous 
events from MOTS (Minitrack Optical Tracking System) and laser systems on 
GEOS-I and II, the BC-4 worldwide camera network on PAGEOS, and local datum 
survey ties and baselines. The reduced condition equations are developed in this 
appendix and a case is considered for the treatment of correlated observations. 


1. Technique for the Normal Equations 

The mathematical analysis leading to the formation of normal equations for 
the geometric adjustment of coordinates of tracking stations is based on the fol- 
lowing type of events: 

1. Two cameras observe the satellite simultaneously 

2. Three cameras observe the satellite simultaneously 

3. Four cameras observe the satellite simultaneously 

4. Two cameras and one laser observe the satellite simultaneously. 

Condition equations resulting from a given set of simultaneous observations 
are of two types: 

• Coplanarity equation, which requires that two observing stations and 
their directions to the satellite lie in the same plane. 

• Length equation, which requires that the satellite position satisfying the 
two-station coplanarity relationship also agrees with the range from a 
third station. 

Corresponding to each event condition equations of the following form are 
used: 


( 1 ) 

i j 

ERECTING page blank NOT FILMED 
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where a. , b., and c are known constants derived in the subsequent sections, 
c is the discrepancy in the condition equation 

v . are the unknown residuals (adjusted minus observed values) 

x. are the unknown corrections to stations' Cartesian coordinates 
(adjusted minus initial values) 

m is the number of observed quantities 

n is the number of unknown coordinates. 

The number and types of condition equations for events 1 through 4 above 
are as follows: 

• For a two-camera event, one coplanarity equation is used. 

• For a three-camera event, three coplanarity equations are used. 

• For a four-camera event, five coplanarity equations are used. 

• For a two-camera, single-laser event, one coplanarity equation and one 

length equation are used. 

The number of condition equations for each event corresponds to the number of 
observations less three, since the observation equations are reduced to a form 
where the three satellite position coordinates are eliminated. Each observing 
camera contributes two observations in an event and an observing laser contrib- 
utes one observation. Each of the coplanarity equations for an event involves 
distinct pairs of observing stations. 

Additional condition (constraint) equations specify coordinate differences 
and also take the form of equation (1). Constraints are treated statistically, 
similar to observation equations, where values and accuracies are obtained from 
a priori information based on datum survey. Two types of constraint equations 
are applied: 

• Distance equations (baselines), which require the distance between two 
stations to remain near a given value 

• Coordinate-shift equations, which require the differences between co- 
ordinates of two nearby stations to remain near a given value. This 
constraint is used to connect the stations in the geometric geodesy with 
those in the dynamic satellite geodesy. 

For a geometric only solution a third type of statistical constraint may be 
applied on individual station coordinates in order to fix the origin of the system. 
A priori values for these constraints should be taken from a different source 
than datum survey such as from a previously determined geocentric solution in 
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a center of mass reference frame. This constraint is not used in the combina- 
tion solution with dynamic satellite geodesy. 

For each event (or constraint) k, denote the associated condition equations 
of the form ( 1 ) in matrix notation as 

AfeVfc + BfeX + Q, = 0, (2) 

for which an example of the dimensions and elements of the matrices are given 
below. Minimizing Q below w.r.t. the unknown station coordinates in X and 
residuals in will lead to the formation of the normal matrix equation. The 
form Q is 


K 

Q = (V iVfc - * V * < 3 > 

where W 4 is the diagonal weight matrix for the observations in and each 
is a column vector of Lagrangian multipliers corresponding to the number 
of condition equations in event k . The resulting normal matrix equation to be 
combined with the gravimetric and dynamic satellite geodesy systems is 


<4) 

where 

K 

(5) 

k= 1 

% = (A^' A .{) (6) 


The largest dimension of is 5 x 5 corresponding to event of type 3 where 
there are 5 coplanarity equations of condition. This case has 8 observations and 
the dimensions of Afc., V&, Bfo, and C& are respectively 5 x 8, 8 x 1, 5 x N x , and 
5 x 1 . N x is the total number of all station coordinates and Bfj (5 x N x ) would 
contain for each of the 5 rows only 6 non-zero elements, corresponding to the 
b. coefficients in (1) for each distinct coplanarity equation involving two observ- 
ing stations. Each row of has 4 non-zero elements corresponding to the 
coefficients in (1) , associated with the two observing stations in each coplanarity 
equation . 
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By employing a suitable set of constraints including those that fix the origin, 
N may be set equal to zero and a geometric only solution for X can be derived 
from (4). 

Condition equations for coplanarity, length, and constraints are developed 
in sections 2 through 5 and section 6 treats a case for correlated observations. 


2. Coordinate System 

Camera observations in a and 8 are transformed from right ascension a 
and declination 8 to earth-fixed angles /3 and y. The conversion of a and 8, 
as corrected for precession, nutation, and polar motion, to the angles fi and y 
is straightforward. The topocentric angle y is measured with respect to the 
equatorial plane and is equivalent to 8 , i.e., y - 8 . The angle /3 is measured 
from the Greenwich meridian in a plane parallel to the equator and is 

J3= a - GHA (7) 

where GHA is the Greenwich Hour Angle at the epoch of the observation. 


3. Coplanarity Equation 

The coplanarity equation requires that the volume of the parallelepiped 
defined by the two station-to-satellite vectors and the station-to-station vector 
and their respective errors be zero. The two station-to-satellite vectors are 
defined in the local terrestrial coordinates as 


where 


p. = u. i + v. j + w. K 

*1 1 t J 1 


( 8 ) 


u. = cos y. cos B. 

i y i ' i 

v. - cos y. sin (i = 1, 2) 


w. = s i n y. 

The station-to-station direction vector p 3 is similarly defined in spherical 
coordinates by use of 





0 < Aj < 277 


(9) 
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7 3 = tan -1 


Z 2 - Z 1 


((* 2 - x ,r + ( y 2 - yi > ) 


77 . .77 

£7, < ~ , 

2 \ 1/2 2 3 2 


where x, y, z are the Cartesian coordinates and the range between the station 


r 3 = (( x 2 - x ,) 2 + ( y 2 - y ^ 2 + ( z 2 - z ,) 2 ) 1/2 

The volume of the parallelepiped defined by these vectors (p t , P 2 , P 3 ) is given 
by their triple scalar product, which is the determinant 


cos 7 . cos /3 cos 7 cos j3 cos 7 cos f3 

1 1 2 2 '3 o 

F Q - cos sin ySj cos y 2 sin ^ cos y 3 sin 


sin 7 


sin y 2 


sm y 


and the adjusted volume through linear expansion is 


F = + AF - 0 


The coefficients of the expansion are then given by 


o r q 

a, = = cos y, sin y 2 cos y $ cos(/S 3 - fS^ - cos y t cos y 2 sin 7 3 cos(/S 2 - ( 13 ) 


2 By, 


- COS 7j cosy 2 cos y 3 sin(/S 3 - /3 2 )- sinyj cosy 2 sin y 3 sin(/3 2 - 


s in 7 j s in y 2 cos y 3 sin(/5 3 - /S ) 


a 3 s TJ~ = cos 7 2 [cos y l sin y 3 cos(/3 2 - - sin y x cos y 3 cos(/3 3 - &,)] 
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= - cos 7 j cos y 2 cos y 3 si n(^ 3 - - sin y 1 sin y 2 cos y 3 sin(/3 3 -Q) 

- cos y i sin y 2 siny 3 sin(/3 2 - 




— = cos y 3 [s in y % cos y 2 cos(/? 3 - /3 2 ) - cos y x sin y 2 cos(/3 3 - /^)] 


( 16 ) 


(17) 


b 


2 


d l± 

^3 


= cos y x cos y 2 cos y 3 sin(/? 2 - /? 1 ) - sin y x cos y 2 sin y 3 sin(/3 3 - /3 0 ) 
+ cos y x sin y 2 siny 3 sin^ 


(18) 


Since , yj , and y 2 are observations, A/3 X , Ay x , A/3 2 > and A y. are residuals 
and are designated v x , v 2 , v 3 , and v 4 , respectively. A/3 3 and Ay 3 are the inter- 
station direction adjustments. The variables to be solved for are corrections to 
the stations Cartesian coordinates. The transformation of unknowns from inter- 
station direction to Cartesian coordinate corrections are given by equation 26. 
Then there results an equation of the form of (1). 


4. Length Equation 

The length equation is developed for two cameras and a laser DME observ- 
ing the satellite simultaneously. Assume the existence of two cameras (A and 
B), the laser DME (L), and the satellite (S), where directions from the cameras 
to the satellite are observed simultaneously (A to S and B to S) and a range is 
observed at the same time from L to S. These quantities and auxiliary vectors 
and angles are shown in Figure A-l. Assumed values of coordinates of the 
cameras and the laser system are used to calculate initial estimates of the 
directions and distances between the cameras and the laser. By taking scalar 
products of the station-to- station and station-to-satellite vectors the cosines 
of the angles, £ , 17 , and £ are obtained as follows: 

cos g = ~P{P 2 - sin y x sin y 2 + cos y x cos y 2 cos(/3 2 - ( 19 ) 

cos rj = ~p{P 3 - s in y x sin y 3 + cos y x cos y 3 cos(/S 3 - (20) 


24 




cos £ = p^p 4 = sin y 2 siny 4 + cos y 2 cos y 4 cos(/S 4 - /3j 


( 21 ) 


where /3 t , y x , fi 2 , and y 2 are directions to the satellite, 

/3 , y 3 are the inter- station angles for the two cameras, and 

/3 , y are the inter-station angles for one camera and the laser 

system. 

From Figure A-l the law of cosines will give, corresponding to the laser 
length s 


F = r 2 + b 2 - 2br cos £ - s 2 = 0, 
and the law of sines will give for b above 


( 22 ) 


, sin ri 

b = a 

sin 

Through the use of (19) through (21) we expand F linearly about the values of 
a, r, /3 3 , y 3 , , y 4 , obtained from the initial station coordinates, and the values of 

» ^2» y 2 ' and s o f rom the observations. Then we have using differentials 
as adjustments (d s A) 


F = F, 


3 F 


3F 


OF 




+ ds = 0 


Divide F through by q = 2(b - r cos £) and denote the result by 


a i d ^l + a 2 dy i + a 3 d ^2 + a 4 d ">2 + a s ds + b l d ^3 + b 2 d:/ 3 + b 3 da + b 4 d ^ 4 + ^d^ 

+ b g dr + c = 0 

where C = F 0 /q, and the differentials on the a i coefficients are the observation 
residuals v 4 for i = 1 to 5. This represents the laser length equation. The co- 
efficients and C are evaluated from the initial values, where F 0 is obtained from 
the misclosure of (22) and the coefficients in (23) are as follows: 
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\ = Pj/sin £ -P s /sim7 


b i = 


P s /sin rj 


( 24 ) 


a 2 ” P 2 / sin £ " P 6 /sin V 


b 2~~ 


- P 8 /sin 77 


b 3 = 


sin 77/sin g 


a 3 = - P 1 /sin £ - P g /sin $ 
a 4 = P 4 /sin^-P 10 /sin £ 
a = - s/(b - r cos £) 


b 4 ~ 


Pg/sin i 


- P 12 /sin l 


b 6 = 


(r - b cos £)/(b - r cos £) 


and where the P's are given as 

Pj = b cot £ [cos y 2 sin(/3 2 - ^)] cos y x 

P 2 = b cot £[cos y x sin y 2 - sin y x cos y 2 cos(/? 2 - ^)]- 


P = - P 

r 3 l 


P 4 = b cot ^[sin 7j cos y 2 - cos y x siny 2 cos(/3 2 - /3 1 ^)3 


P e = 


a cos 77 


5 sin£ 


[cos y t cos y 3 sin(/? 3 - Z^)] • 


p = ~- t ?° S 1 [cos y. sin 7 , - sin 7 cos 7 cos(/3 - /3 )] 

6 sin£ 13 13 3 1 


P 7 = 


-P* 


P 8 = 


a cos 77 
sin 


[sin 7 vCOS 7^ - cos 7j s in 7 3 cos(/S 3 - Z^)] 


P 9 = 


= ( br sin L\ [ cos y COS 7 S in(/3 - 8)] 

\b - r cos i/ 1 4 4 1 
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p 


p 


p 


10 


/ br sin £ \ 
\b - r cos 


n = - p * 


I br sin £ \ 
\b - r cos £/ 


[cos y 2 sin y 4 - sin y 2 cos y 4 cos(/3 4 - 0 2 )] 


[sin y 2 cos y 4 - cos y 2 sin y 4 cos(/3 4 - f3 2 )] 


In order to obtain the desired form (1) for Cartesian station coordinates, the 
coplanarity and length equations are transformed from 7,/?, r variables to x, 
y, z variables by using the relationships 


x 2 - Xj - r cos 7 cos /S 
y 2 - y 1 = r cos 7 sin /3 

z 2~ z l = r sin 7 

Differentiating these expressions yields 


dx 2 - dx t 


-r cosy sin/? -rsinycos/3 cos y cos /3 


H “1 

d/3 

dy 2 - dy t 

= 

r cos y cos /3 -rsinysi n j3 cosy sin/3 


dy 

dz 2 - dz t 


0 r cos y sin y 


dr 


Inverting Equation 26 produces the transformation 


d/3 


-sin jB 

cos jB 

0 


dx 2 - dx t 


r cosy 

r cos y 




dy 

= 

-sin y cos j3 
r 

-sin y s in /3 
r 

cos 7 
r 


dy 2 - dy } 

dr 


cos y cos jB 

cos y sin ft 

sin y 


dz - dz 

L 2 J 


(25) 


(26) 


(27) 
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5. Condition Equations for Constraints 

The three types of condition equations are: (1) coordinate equations, 

(2) distance equations, and (3) coordinate- shift equations. As indicated pre- 
viously constraints for station coordinates, distances (baselines), and coordinate 
shifts (local datum ties) are based upon a priori information. Such information 
is treated statistically as in the case of satellite observations, and hence weights 
are applied corresponding to a priori errors. 


5.1 Coordinate Equation 

Assume the input coordinates of the ith station are coordinates for which 
a priori information is available. Let the input values of X. , Y . , Z . be X. o , 

Y , Z. , and denote the adjusted coordinates as 

dx. = Xl - X io - x._ 2 

d Vi=Y. -Y io *Vi < 28) 

dz i = Z i - Z io s x i 


then the constraint equations in the form of equation (1) are simply 


V. A 

- X. 

= o 

i-2 

i-2 

V. „ 

- X. , 

= 0 

i-2 

1-1 


V. 

L 

- X. 
1 

= 0 


(29) 


5.2 Distance Equations (Baselines) 

The condition equation for the baseline distance q between the rth and sth 
ground stations is 

-dq + cos y cos /3(dx s - dx r ) + cos y sin J3( dy s - dy r ) + siny(dz s - dz f ) = 0 (30) 

where y and /3 are used as in (9) and (10) and the differentials are the unknown 
station adjustments. The adjustment 

dq = q- q c = (q~<^ ) ) - (q^ - q Q ) 

= v - C (31) 
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where q is the solution value, q c is the computed value based upon the initial 
station coordinates, and q 0 is the a priori value for the constraint. 

In conformity with the previous notation the terms dx s , dy s , dz s , dx r , dy r , 
and dz r are replaced by x s _ 2 , x s _j, x s , x r _ 2 , x r _ t , s r , and dq by (31) to obtain 

-v + cos yc os/3(x s _ 2 - x r _' 2 ) + cos y sin/3(x s _ 1 - x^) + sin y(x s - x p ) +C = 0 (32) 


5.3 Coordinate-Shift Equations 

For two nearby stations denote the difference in coordinates as 



(33) 


Since the results are similar for the three equations we will treat just one equa- 
tion of condition. For the x component the differential is 


- dD x + dx 2 - dXj = 0, 

and as in the case of the distance equation (31) 


(34) 


dD x = 


v -C 




(35) 


where D x is the unknown difference, D X() is obtained from local survey station 
coordinates, and D Xc is computed from the initial input of the station coordinates. 
For the ith and jth stations, using the notation of the form (1) where differentials 
are replaced by corrections x q and x p , the condition equations for (33) are. 


- V 

+.X. ^ 

- X. ^ 

+ C 

* 

J~2 

*-2 

X 

- V- 

,+ X. - 

— X. „ 

+ C 

y 

rl 

l - 1 

y 

- V 

z 

V X. . 

v * j / 

- X. 
1 

+ C 2 


0 

0 

0 


(36) 
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6. Correlated Observations 


The model described above was developed for camera systems that observed 
simultaneously the flashing lamps on GEOS-I and II, and then the model was 
employed to include the BC-4 camera network that observed the PAGEOS satel- 
lite. A BC-4 photograph taken on PAGEOS by an observing station, s, was re- 
duced to 7 time points (It = 1 to 7) of satellite observation angles The 

reduced observations y|, and /3| are correlated separately in each type among 
the points It = 1 to 7. The modeling for the correlated observations is presented. 

Consider 7 events of the type (1), (2), or (3) described in section 1, where 
respectively 2, 3, or 4 stations (S - 2, 3, or 4) observe the satellite simultaneously 
at each of the 7 reduced photographic points It = 1 to 7. Thus for each event 
there are 2S simultaneous observations, namely (7^,/^) for s = 1 to S. Let p 
denote this configuration of S stations and 7 events, then for each p there are 
7 sets of matrix condition equations of the form (2). Denote these as 

A v +B X + C - 0 (37) 

P P P P 

where by row partitioning for It = 1 to 7 

v [v & 

C p =[c|] (38) 

B p = [Upl- 
and A| lies along the diagonal submatrix path of A p (with zero submatrices 
for off diagonal blocks) 


A p =[DlAGA£] (39) 

The submatrices v£, B|, and A£ are given as before in (2) for a particular 
event, but here the event type for S = 2, 3, or 4 stations is fixed for the 7 events 
for a given configuration p. 

Denote the variance-covariance matrix of the observation errors as 

P n - E(Y vj) (40) 

p v p p ' 

where V p corresponds to the observation errors (noise) in V p . 
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It is of interest to compare the matrix M p derived from 7 events to that of 
given in (6) for a single event k. Take the case of S = 4 for which the dimen- 
sions Aft,, Vf^ , and were given under (6) respectively as 5 x 8, 8 xi, and 
8x8, and for which there were 8 observations and 5 coplanarity equations of 
condition from the 4 observing stations in a given event. Consider 7 events of 
the same type as in the configuration p for S = 4, denote as M|. for k = 1 
to 7, and assume correlations are absent as in the previous modeling. Then 
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0 


H» = 


M, F 


MP 




■ CDIAGM|] 


(45) 


35* 35 


where 


4 = A£(«1> M (A|>\ 

b A 5 

and since correlations are assumed absent here 


(46) 


P p = [DIAG(W|) _1 ], 

With correlations the same diagonal blocks in (45) arise for M in (44) since 
in each event h all observations are uncorrelated, but similar off diagonal blocks 
also exist which is now shown. Using the submatrices VjP for V p in (38) and 
dropping the superscript p on the submatrices, then the variance-covariance 
matrix in (40) becomes 

P p = E(V p V p ) = CE(%vJ)] for k = 1 to 7 and l = 1 to 7, ( 47 ) 

which corresponds to 49 sublocks or submatrices in P p . For a given event k 
and i the only covariances occur when the station s and the angle or 
are the same. Denote the observation errors for a given k as (where T = true, 
o = observed) 


^(T) - yj>( 0) 

- ?i<P> 


V 


r|(T> - y|(0) 
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s [v.^] , i = 1 to 2S> 


( 48 ) 



then for a given k and £ 


E <W= ^ik v ii )] 2 s« 2 s j = 1 to 2S , 

and from the definition of the correlations 

E(v. £) = o-?. (k, ■i ) = 0 for i ^ j for all k andZ 

= 0 - 2 . (k, -i) for covariances h \ ^ , 

= 0 - 2 . (&, k) for var iances h = 't. 

Thus, in each sublock of a given k and 4, the off diagonal elements 


(49) 


E(V^vJ) = (pi AG (A, 1 % 


2SX 2S 


(50) 

are zero and 

(51) 


s D. 


U 


and 


P p = [D^i (k l = 1 to 7), 


(52) 


Denote 


WT — TlUo ol ■ (k. 




and with use of (38), (39), and (52) in (44) then by (53) 

~ (2 s- 3) x (2 S~ 3) 

Now 

\k = Wk A l 

'=W A l = % 


(54) 

(55) 


which is the same as in (46) for the uncorrelated case. 

The block form [M^] for M , where M.jy- is given in (54), provides a con- 
venient method for the computations of M . However in the present case of cor- 
related observations there are 49 such blocks, whereas only the 7 diagonal blocks 
are computed for uncorrelated observations as in (45) or (55). The largest 
inverse matrix M p to be inverted occurs for the case of S = 4 and which has 
dimension 35 x 35, whereas previously for the case of uncorrelated observations 
the largest matrix was 5x5. Correlations are generally large among the 
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reduced observations of a photograph. Thus the geometric normal equations 
should be analyzed further to investigate the overall effect of the correlation on 
the final combination solution. 
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